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1. SUMMARY 

This report describes a numerical method for generating three- 
dimensional grids for general configurations. The basic method involves 
the solution of a set of quasi-linear elliptic partial differential equations via 
pointwise relaxation with a local relaxation factor. It allows specification of 
the grid spacing off the boundary surfaces and the grid orthogonality at the 
boundary surfaces. It includes adaptive mechanisms to improve smoothness, 
orthogonality, and flow resolution in the grid interior. 

2. INTRODUCTION 

Three-dimensional computational fluid dynamics codes require compu- 
tational grids with suitable resolution, smoothness, and orthogonality. High 
grid resolution allows complex flow physics to be modelled near shocks and 
in shear layers. Smoothness of the metric data prevents the flow solution 
from being dominated by truncation error in the metric coefficients. Grid 
orthogonality at the boundaries simplifies and improves the accuracy of any 
boundary condition involving normal gradients. 

This report describes a numerical method for generating three- 
dimensional grids for computational fluid dynamics codes with the potential 
and the ability to satisfy the foregoing requirements. The basic method is 
general and involves the solution of a quasi- linear elliptic partial differ- 
ential equation via pointwise successive relaxation with a local relaxation 
factor. The governing equation contains forcing functions that depend upon 
the boundary point distribution and the boundary surface gradient. The 
method allows specification of the grid point distribution on the boundary 
surfaces, the grid spacing off the boundary surfaces, and the grid orthogo- 
nality at the boundary surfaces. It includes adaptive mechanisms to improve 
smoothness, orthogonality, and flow resolution in the grid interior. 
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The mechanism of adaptation contains two branches: cumulative move- 
ment and non-cumulative movement each with its own characteristic. 

In non-cumulative action, the amount of adjustment made on the nodal 
point distribution is governed by a penalty type of behavior, with fixed upper 
and lower limits, on the forcing functions. 

The system appears to be very stable numerically. The overall accuracy 
of the system, in terms of residues of the physical aspects under adaptation, 
is limited by the strength of the penalty function. 

In cumulative action, the adjustments are guided by physical con- 
straints. Movement of the grid will not stop unless a homogeneous condition 
is reached. There are no directly predetermined upper or lower limits on 
the forcing functions with this technique. 

The strength of the adaptive mechanism is determined by a scaling 
factor, which in turn controls the magnification of disturbances at each 
iteration. The behavior of this technique resembles an explicit scheme, the 
scaling factor should be small enough to assure stability of the system, and 
yet large enough to achieve better accuracy. 

This paper includes a discussion of the mathematical formulation of the 
method and some representative results. 

3. MATHEMATICAL FORMULATION 

The quasi-linear elliptic governing equation is taken from reference 1 as 
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The metric tensor components g' 3 and g kk in equation (I) are defined 


as 
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The P* in equation (1) are the forcing functions specified by the user. 
They represent one-dimensional stretching in each index coordinate direc- 
tion. 

Equation (1) can be rewritten using matrix notation as 

A + B P = 0 (2) 

where 
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Equation (2) can be solved for P* on the boundaries as 

P 0 = -Bo 1 A 0 (3) 

where the subscript “0” indicates values on the boundary. Tangential deriva- 
tives for terms on the right hand side of equations (3) are determined by 
applying standard difference formulas to the prescribed boundary point dis- 
tribution on the surface. Normal derivatives are determined by specifying 
the first normal derivative equal to the desired spacing off the boundary (for 
detail, see reference 2) and using the approximation 

<3 2 zq _ 2(xi - xq) 2_ dip 

3n 2 (An) 2 Arc drc 

where n indicates the normal direction, the subscript “0” indicates values on 
the boundary, and the subscript “1” indicates values one point away from 
the boundary. 

Once the boundary values are known, the interior values of Pjt can be 
determined using 
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The value of Po t , , represents the fc-th component of the P 0 vector on 
the minimum boundary surface in the /-th direction. The value of Po fc , i 
represents the fc-th component of the P 0 vector on the maximum bound- 
ary surface in the /-th direction. The functions a , /5, and 7 have subscript 
notation similar to that of Po- The a function represents polynomial ex- 
trapolation from a controlled boundary using a constant factor C a . The 7 
function represents the combined effect of the /3 functions, which represent 
power-law factorization with constant exponent C$ to control the depth of 
influence away from a controlled boundary. 

The values of the forcing functions can be modified for improved 
smoothness by using 


P'k = QPk 


where 


6 = 1 - tanh(C«,(l - a c> »)) 


* = J/(\\ '1 II • II h || • || h ||). 

The constants C e , and Cg 7 define the rate and order of the adaptation. 
The variable a is a measure of the shear of the grid, with J representing the 
Jacobian and / 1 , / 2 , I 3 representing the grid cell lengths in each direction. 

A measure of the local orthogonality of the grid can be defined as 
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The constant C$ defines the order of the adaptation. A one-dimensional 
variational form of equation (4) can then be written as 
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If a flow error index E is computed by the flow solver such that E > 0, 
a measure of the local flow resolution can be defined as 

* = J{1 +f(E)). (5) 

where / represents a suitable functional of E. 

A one-dimensional variational form of equation (5) can then be written 
as 
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The values of the forcing functions can be modified for improved local 
orthogonality and flow resolution using 


P' k = {0 + • 6 x )Pk Fk-Pk> 0 

( 6 ) 
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where 9x, similar to 9, is the skewness index for the adaptive constributions, 
and 

A fc = C Al tanh(C A5 |F fc ! c ^) 

F k = {W^' k +W^' k )/m + W^) 

with the constants Cx x , C>, 2 , and determining the range, rate, and 
order of the adaptation for the non-cumulative mechanism. The variable F 
represents a weighted combination of the skewness and flow error variations. 

When cumulative adjustment is desired, the focus must be on the dis- 
turbances of F k instead of F k itself. Thus let 

A F k = F { k n ~ 1] - F { k n) 
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X k ‘ = tanh(C> 3 A£) when AjJ. > 0 

A“ = 0 when X k < 0 

and then 


Afc = A*. n_1) + X’ k (7) 

Equation (7) is used in conjuction with equation (6) for evaluating the 
proper value of the forcing functions at the N t /, step. 

4. RESULTS AND DISCUSSION 

Two groups of test cases using a turbine blade geometry are attached. 
Each of these cases represents a particular capability of the technique. 
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Within each group, the nodal point distribution on the boundary lines and 
the basic values of the forcing functions are constructed in exactly the same 
fashion. 

Figure 1 illustrates the effect of adaptation according the the skewness 
of the grid only. Highly sheared grids (around the trailing edge) are suffi- 
ciently improved compared to results without this adaptation (not shown). 

Figure 2 shows the effect of the non-cumulative adaptation according 
to the skewness and the flow property with constant error index everywhere 
(equivalent to requiring a zero variation in the Jacobian). The result has 
indicated that the nearly collapsed grid lines (from Figure l) at the inlet 
corner are greatly improved whereas the grid elsewhere is well preserved. 

Figure 3 illustrates the effect of the non-cumulative adaptation accord- 
ing to the skewness and the general orthogonality of the grids. The result 
shows the tendency of the expansion of a good grid and the suppression of 
a sheared grid. 

Two representatives of the three-dimensional capability of the current 
technique are illustrated in Figure 4 and 5. Figure 4 shows a C-grid configu- 
ration with adaptation according to the skewness of the grid only. Figure 5 
shows an H-grid configuration with the adaptation according to the skewness 
and the non-cumulative orthogonality of the grid. 

Figure 6, 7, and 8 illustrate the effects of the solution adaptive mecha- 
nisms. The error index used here is a fixed artivicial pattern with 

E — ((J - l)/JHalf)* when J< JHalf + 1 

= (( Jmax — J)/ J Hal})* when J > JHalf + 1 

where J Hal f = ( Jmax — l)/2 + 1. 

The error functional / used is 

1 + f{E) = e AE 

with a constant A. 

Figure 6 represents the optional grid with a skewness adaptation only. 
Figure 7 shows the effect of the non-cumulative solution adaptive mechanism 
with A = 10. Figure 8 illustrates the effect of the cumulative solution adap- 
tive technique with A = 2. The results have indicated a positive response on 
the redistribution of the grids. The adjustments achieved by the cumulative 
mechanism are more profound than that indicated by the non-cumulative 
mechanism. 

The stability of the adaptive technique is satisfactory. Failure to reach 
convergence is possible, however, particularly when the spacing requirements 
of the grids are extremely small and when the physical configuration forces 
the appearances of highly sheared grids in the interior. These difficulties 
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can be relaxed, only to a degree, by the specification of a suitable strength 
for the adaptive schemes. 

The solution adaptive technique introduced here is still at its early 
developing stage. The adaptation according to a real slow property has not 
yet been tested. Although current test results have shown positive response 
on the grid movements to the adaptive mechanisms, they also appear to 
indicate a major weakness of this technique, that is, the adjustments on the 
grid are ultimately limited by the basic requirements, which are specified 
by the user. The residue of the property under adaptation can only reach 
a relative minimum, not an absolute minimum, in order to maintain the 
numerical stability. 
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Fig- 4 Three dimensional C-grid configuration. 



Fig. 5 Three dimensional H-grid configuration. 
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